# ***************************************************************************
# *                                                                         *
# *   Copyright (c) 2016 - Qingfeng Xia <qingfeng.xia iesensor.com>         *
# *                                                                         *
# *   This program is free software; you can redistribute it and/or modify  *
# *   it under the terms of the GNU Lesser General Public License (LGPL)    *
# *   as published by the Free Software Foundation; either version 2 of     *
# *   the License, or (at your option) any later version.                   *
# *   for detail see the LICENCE text file.                                 *
# *                                                                         *
# *   This program is distributed in the hope that it will be useful,       *
# *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
# *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
# *   GNU Library General Public License for more details.                  *
# *                                                                         *
# *   You should have received a copy of the GNU Library General Public     *
# *   License along with this program; if not, write to the Free Software   *
# *   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  *
# *   USA                                                                   *
# *                                                                         *
# ***************************************************************************

from __future__ import print_function, absolute_import

import multiprocessing  # used in BasicBuilder to detect CPU cores
from .utility import *
from .utility import _debug

"""
incompressible flow with RAS and LES turbulence setup,
transient, dynamic meshing support should be implemented in this class
boundary conditions naming modelled after Ansys Fluent/CFX
Boundary type and subtype must consistent with definitioin in FemConstraintFluidBoundary.cpp of FreeCAD FemWorkbench
"""
supported_boundary_types = set([
'wall',  # default
'interface',
'inlet',
'outlet',
'freestream',
])

# wall is the default boundary type
supported_wall_types = set([
'fixed', # fixedValue/noSlip
'slip', # perfect 100 percent slip wall
'partialSlip', #
'moving', # movingWallVelocity
'rough', # for turbulence flow and using special wall function
#'translatingWallVelocity' #
#'atmBoundary', # modelled atomspheric boundary layer on earth surface
])

"""
constraint interface boundary:   `ls $FOAM_SRC/finiteVolume/fields/fvPatchFields/constraint`
No need to pair/coupled with other boundaries
"""
_constraint_interface_types = set([
'symmetryPlane',  # planar symmetric, zeroGradient for all
#'symmetry',  # for any (non-planar) patch which uses the symmetry plane (slip) condition
'empty',   # 2D domain, default name is "frontAndBack", axis of axis symmetric
'wedge',   # wedge front and back for an axi-symmetric geometry
])

"""
coupled interface boundary: need to be paired with another boundary
`ls $FOAM_SRC/finiteVolume/fields/fvPatchFields/`
"""
_coupled_interface_types = set([
'cyclic',  # Enables two patches to be treated as if they are physically connected, pressure boundary for geometry in grid pattern
#'cyclicAMI', #like cyclic, but for 2 patches whose faces are non matching; used for sliding interface in rotating geometry cases.
#'processor',  # inter-processor mesh decomposition boundary, automatically generated by domain decompose tool
'coupled', # exchange boundary value with external solver by file
#'mapped'
])
supported_interface_types = _constraint_interface_types | _coupled_interface_types

"""
User defined function for boundary condition like non-uniform velocity profile is not supported
CSV table,  codedStream, timeVarying, variable, could be used in manually case editing
see <http://cfd.direct/openfoam/user-guide/boundaries/>
"""
supported_inlet_types = set([
'totalPressure',
'staticPressure',
'massFlowRate', # flowRateInletVelocity
'volumetricFlowRate', # flowRateInletVelocity, flowRateInletVelocity
'uniformVelocity', # it is called "fixedValue" in OpenFOAM dict file
])

supported_outlet_types = set([
'staticPressure',
'uniformVelocity',
'outFlow'  # it is called 'inletOutlet' in OpenFOAM
])


def getDefaultSolverSettings():
    """ default sovler settings with the first letter as lower key
    corresponding to CfdSolverFoam object' s properties
    """
    return {
            'parallel': False,
            'compressible': False,
            'nonNewtonian': False,
            'transonic': False,
            'porous':False,
            'dynamicMeshing':False,
            'buoyant': False, # body force, like gravity, needs a dict file `constant/g`
            'gravity': (0, -9.81, 0),
            'transient':False,
            'turbulenceModel': 'laminar',
            # CSIP team contributed feature,
            'potentialInit': False, #  new property inserted into CfdSolverFoam
            #
            'templateCasePath': None,  # New feature in 2020, control how case is created
            'caseCreationMode': "fromScratch",
            # heat transfer specific properties
            'heatTransfering':False,
            'conjugate': False, # conjugate heat transfer (CHT)
            'radiationModel': 'noRadiation',
            #'conbustionModel': 'noConbustion',
            #
            #'multiPhaseModel': 'singlePhase' # twoPhase, multiphase
            #'missible': False, # two species can be mixed perfectly if missible == True
            }  # containing all setting properties

def _getMultiphaseSolverName(settings):
    # solver list: https://github.com/OpenFOAM/OpenFOAM-5.x/tree/master/applications/solvers/multiphase
    # interFoam: 2 incompressible, isothermal immiscible fluids using a VOF (volume of fluid)
    # multiphaseModel: mapping to controldict of 
    # interfaceTrackingMethod for different phases: VOF, Euler
    # spcie count: 2, 3, n
    # special, cavitity, film, phaseChange, twoLiquidMixingFoam
    if 'compressible' in settings and settings['compressible']:
        if settings['multiphaseModel'] == 'VOF':
            #interMixingFoam  Solver for 3 incompressible fluids, two of which are miscible, using a VOF
            pass  # todo

def _getSolverName(settings):
    """ user may also provide the solver exmplicitely
    if case is builted from solver, solver should be extracted from ControlDict['solver']
    """
    if settings['dynamicMeshing']:
        assert settings['transient']

    # LES model relies on the case template for wall functions and turbulenceProperties
    if settings['turbulenceModel'] in LES_turbulence_models:
        if settings['compressible']:
            return 'rhoPimpleFoam'
        else:
            return 'pisoFoam'
    if settings['turbulenceModel'] == "DNS":
        return 'dnsFoam'
    if settings['turbulenceModel'] == "invisid":
        return 'potentialFoam'
    #
    if 'heatTransfering' in settings and settings['heatTransfering']:
        # for highly compressible flow, just use compressible flow solver which must solve temperature filed
        if settings['dynamicMeshing'] or settings['porous'] or settings['nonNewtonian']:
            print("Error: no matched solver for heat transfering, please develop such a solver")
            raise NotImplementedError()
        if settings['transient']:
            if settings['conjugate']:
                return 'chtMultiRegionFoam'
            if settings['compressible']:
                return 'buoyantPimpleFoam'  # natural convection of compressible
            else:
                return 'buoyantBoussinesqPimpleFoam'  # natural convection of incompressible
        else:
            if settings['conjugate']:
                return 'chtMultiRegionSimpleFoam'
            if settings['compressible']:
                return 'buoyantSimpleFoam'
            else:
                return 'buoyantBoussinesqSimpleFoam'
    #
    if 'multiphaseModel' in settings and settings['multiphaseModel']:
        # this should be a python dict
        print("Error: multiphase model case builder is not implemented yet in _getMultiphaseSolverName()")
        raise NotImplementedError()

    if 'compressible' in settings and settings['compressible']:
        # selection of compressible flow solver also depends on Ma range,  liquid compressible is a special case
        # Incomressibile (Ma < 0.3), high Sub-sonic ( 0.3 < Ma < 0.7 ), Trans-sonic ( 0.8 < Ma < 1.0) and super-sonic ( Ma > 1.0), hypersonic (Ma > 5)
        if "compressibilityRegime" in settings:
            print('compressibilityRegime: ', settings['compressibilityRegime'])

        if "transonic" in settings and settings["transonic"]:
            if settings['transient']:
                if settings['dynamicMeshing']:
                    return 'rhoCentralDymFoam'   #  Ma>0.7
                return 'rhoCentralFoam'  # this only density-based solver to deal with 
                #  SonicFoam: supersonic and transonic solver using PISO pressure-based, it is transient only solver
                # actually sonicFoam and rhoCentralFoam can deal with "transonic", set "transonic" in fvSolution simple/pimple algo
            else:
                print("Error: steady state transonic solver is not available")
                raise NotImplementedError()
        else:  #subsonic regimes
            if settings['transient']:
                if settings['dynamicMeshing']:
                    return 'rhoPimpleDyMFoam'   # for low Ma: 0.3<Ma<0.7, pressure-based solver
                return 'rhoPimpleFoam'
            else:
                if settings['porous']:
                    return 'rhoPorousSimpleFoam'
                return 'rhoSimpleFoam'
    else:  # incompressible:
        if settings['transient']:
            if settings['dynamicMeshing']:
                return 'pimpleDyMFoam'
            else:
                if not settings['nonNewtonian']:
                    return 'pimpleFoam'  # pisoFoam
                else:
                    return  'nonNewtonianIcoFoam'
        else:
            if not settings['nonNewtonian']:
                if settings['porous']:
                    return 'porousSimpleFoam'
                else:
                    return 'simpleFoam'
    # finally, reached here
    print("Error: no matched solver, print solver setting and raise exception")
    print(settings)
    raise Exception('No solver is deducted from solver settings due to conflicted settings or no such solver implemented')

_VARIABLE_NAMES = {'U': "Velocity", "p": "Pressure",
"T": "Temperature", 
'alphat': "Thermal turbulence",
'alphas': "phase fraction",
'k': "Turbulence Intensity(k)", 'omega': "Turbulence omega", 'epsilon': "Turbulence epsilon",
'p_rgh': "Pressure-rgh"}

#no need for this
_DEFAULT_VARIABLE_BOUNDARY = {'U': "Velocity", 
"p": "Pressure",
"T": "Temperature", 
'alphat': "Thermal turbulence",
'alphas': "phase fraction",
'k': "Turbulence Intensity(k)", 
'omega': "Turbulence omega", 
'epsilon': "Turbulence epsilon",
'p_rgh': "Pressure-rgh"}

def getDefaultBoundarySettings(variable_list):
    pass

def getVariableList(solverSettings):
        # density variable/field 'rho' will be automatically created by the compressible solver if not present
        vars = ['p', 'U'] + getTurbulenceVariables(solverSettings)
        if solverSettings['buoyant']:
            vars.append("p_rgh")  # commpressible flow?
        if solverSettings['dynamicMeshing']:
            vars.append("pointDisplacement")  # only 6DoF motion solver needs this variable
        #multiphase flow is not yet supported
        return list(set(vars + _getThermalVariables(solverSettings)))

def _getThermalVariables(solverSettings):
    """ incompressible air flow must calc var 'T' and 'alphat' for turbulence flow
    heatTransfer should consider buoyant effect, adding var 'p_rgh', even for buoyantBoussinesqSimpleFoam
    radiationModel: G for P1 and fvDOM, IDefault for fvDOM, 
    alphat = turbulence->nut()/Prt; // turbulent Prantl number is in general randomly chosen between 0.75 and 0.95
    """
    if 'radiationModel' in solverSettings:
        radiationModel = solverSettings['radiationModel']
    else:
        radiationModel = 'noRadiation'
    vars = []
    if solverSettings['heatTransfering']:
        if radiationModel == "noRadiation":
            vars += ['T', 'p_rgh'] # p_rgh for buoyant solver far field pressure field
        else:
            print("Error: radiation is not implemented yet")
            raise NotImplementedError()
    else:
        if solverSettings['compressible']:
            vars += ['T']
    if solverSettings['turbulenceModel'] not in set(['DNS', 'laminar', 'invisid']):
        vars += ['alphat']
    return vars


"""Base on tutorials of OpenFoam 4.0, but should works on OpenFOAM 3.0.x to 6.x
not all solver templates are tested
search turbulence model by: `find $FOAM_TUTORIALS -name turbulenceProperties | xargs grep -l buoyantKEpsilon`
"""
_RAS_solver_templates = {
                    'simpleFoam': "tutorials/incompressible/simpleFoam/pipeCyclic",
                    'icoFoam': "tutorials/incompressible/icoFoam/elbow",
                    'pisoFoam': "tutorials/incompressible/pisoFoam/ras/cavity",
                    'pimpleFoam': "tutorials/incompressible/pimpleFoam/TJunction",
                    'pimpleDyMFoam': "tutorials/incompressible/pimpleDyMFoam/movingCone",
                    'porousSimpleFoam': "tutorials/incompressible/porousSimpleFoam/angledDuctImplicit",
                    'nonNewtonianIcoFoam': "tutorials/incompressible/nonNewtonianIcoFoam/offsetCylinder",
                    # compressible
                    'rhoSimpleFoam': "tutorials/compressible/rhoSimpleFoam/squareBend",
                    'rhoPimpleFoam': "tutorials/compressible/rhoPimpleFoam/ras/angledDuct",
                    'rhoPimpleDyMFoam': "compressible/rhoPimpleDyMFoam/annularThermalMixer",
                    # HeatTransferring
                    'buoyantSimpleFoam': "tutorials/heatTransfer/buoyantSimpleFoam/hotRadiationRoom/",
                    'buoyantBoussinesqSimpleFoam': "tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/",
                    'buoyantPimpleFoam': "tutorials/heatTransfer/buoyantPimpleFoam/hotRoom",
                    'buoyantBoussinesqPimpleFoam': "tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling",
                    }

_supported_algorithms = ["SIMPLE", "PISO", "PIMPLE"]

_build_levels = ['scratch', 'template', 'previousCase', 'updatingMesh']

class BasicBuilder(object):
    """ This class contains boundary condition setup for incompressible flow
    transient boundary condition setup is not supported, search openfoam manual: Time-varying boundary conditions
    """
    def __init__(self,  casePath,
                        #meshPath,
                        solverSettings=getDefaultSolverSettings(),
                        #templatePath=None, #detected from solverName which is deduced from sovlerSettings
                        fluidProperties = {'name':'water', "compressible":False, 'kinematicViscosity':1e6, "density": 1e3},
                        turbulenceProperties = {'name':'laminar'},
                        boundarySettings = [],
                        internalFields = {},
                        transientSettings = {"startTime":0, "endTime":1000, "timeStep":1, "writeInterval":100}, # simpleFoam
                        paralleSettings = {'method':"simple", "numberOfSubdomains":multiprocessing.cpu_count()}, # NOTE: check typo should be paralle'l'Settings NOTE: Should we not use scotch for default
                ):
        if casePath[0] == "~": casePath = os.path.expanduser(casePath)
        self._casePath = os.path.abspath(casePath)
        #self._meshPath = meshPath

        self._solverSettings = solverSettings
        self._solverName = self.getSolverName()
        if "templateCasePath" in solverSettings:
            #shoud test existence first, but that does not work for WSL
            self._templatePath = solverSettings["templateCasePath"]
        else:
            self._templatePath = None # self.getFoamTemplate()
        self._solverCreatedVariables = self.getSolverCreatedVariables()

        self._turbulenceProperties = turbulenceProperties
        self._fluidProperties = fluidProperties  # incompressible only
        self._paralleSettings = paralleSettings

        self._boundarySettings = boundarySettings
        self._internalFields = internalFields
        self._transientSettings = transientSettings

    def createCase(self):
        # TODO:        self._solverSettings["caseCreationMode"]
        if self._templatePath:
            createCaseFromTemplate(self._casePath, self._templatePath)
        else:
            createCaseFromScratch(self._casePath, self._solverName)
        self._createInitVarables()  # some high build level leave existing variable files there, just updating boundary values
        createRunScript(self._casePath, self._solverSettings['potentialInit'], self._solverSettings['parallel'], self._solverName, self._paralleSettings['numberOfSubdomains']) # Specify init_potential (defaults to true)

    def _createSolverSolution(self):
        # PyFoam may run into trouble given fvSolution file has C preprocessor command
        # createCaseFromTemplate() should already copy these files
        if self._solverSettings["caseCreationMode"] == "fromScratch":
            createRawFoamFile(self._casePath, "system", "fvSolution", getFvSolutionTemplate())
        else:
            sf = ParsedParameterFile(self._templateCase + "/system/fvSolution")
            of = ParsedParameterFile(self._casePath + "/system/fvSolution")
            of = sf  # todo: unit test
            of.writeFile()

    def _createSolverSchemes(self):
        # createCaseFromTemplate() should already copy these files
        """
        # add new variable (for diff turbulence model) schemes if abscent in source template
        divSchemes
        {
            default         none;
            div(phi,U)      bounded Gauss upwind;
            div(phi,T)      bounded Gauss upwind;
            div(phi,k)      bounded Gauss upwind;
            div(phi,epsilon) bounded Gauss upwind;
            div((nuEff*dev2(T(grad(U))))) Gauss linear;
        }
        """
        if self._solverSettings["caseCreationMode"] == "fromScratch":
            createRawFoamFile(self._casePath, "system", "fvSchemes", getFvSchemesTemplate())
        else:
            sf = ParsedParameterFile(self._templateCase + "/system/fvSchemes")
            of = ParsedParameterFile(self._casePath + "/system/fvSchemes")
            of = sf  # todo: unit test
            of.writeFile()
    
    # TODO: setupCase()  or setup() could be a better name
    def build(self):
        # if case is built from clone/template, this function should not be called, or called with diff build_level
        self.setupBoundaryConditions()
        self.setupInternalFields()
        # build_level 
        self.setupFluidProperties()  # materials properties
        self.setupTurbulenceProperties()

        if self._solverSettings['transient']:
            self.setupSolverControl()
        if self._solverSettings['parallel']:
            # see: http://cfd.direct/openfoam/user-guide/running-applications-parallel/
            self.setupParallelSettings()
            # it is the CfdRunnalbe to mpirun and recompose the result and show result
        if self._solverSettings['buoyant']:
            self.setupGravityProperties()
        if self._solverSettings['dynamicMeshing']:
            self.setupDynamicMeshingProperties()
        self.setupSolverControl() # residual, relaxfactor refValue, refCell etc

        # Move mesh files, after being edited, to polyMesh.org
        #movePolyMesh(self._casePath)  # make trouble in WSL for ln -s command in Allrun script

    def setupMesh(self, mesh_path, scale):
        # create mesh by conversion from other mesh file format
        if os.path.exists(mesh_path):
            convertMesh(self._casePath, mesh_path, scale)

    def updateMesh(self, updated_mesh_path, scale):
        #runFoamCommand('foamCleanPolyMesh -case {}'.format(self._casePath)) #foamCleanPolyMesh v4.0+?
        casePath = self._casePath
        shutil.rmtree(casePath + os.path.sep + "constant" + os.path.sep + 'polyMesh') #  rm -rf dir
        self.setupMesh(updated_mesh_path, scale)

    def check(self):
        """ heatTransfer must has Buoyang = True, and existent of dict
        """
        case = self._casePath
        if self._solverSettings['dynamicMeshing']:
            if not os.path.exists(case + os.path.sep + 'constant/dynamicMeshDict'):
                return "Error: 'constant/dynamicMeshDict' is not existent while dynamcMeshing opiton is selected"
            if not os.path.exists(case + os.path.sep + '0/pointDisplacement'):
                return "Error: '0/pointDisplacement' is not existent while dynamcMeshing opiton is selected"
        if self._solverSettings['porous']:
            return 'Error: Porous flow need temperatuare field, use ThermalBuilder'
        if self._solverSettings['parallel']:
            if not os.path.exists(case + os.path.sep + 'constant/decomposeParDict'):
                return "Warning: 'constant/decomposeParDict' is not existent while parellel opiton is selected"
        if self._solverSettings['buoyant']:
            if not os.path.exists(case + os.path.sep + 'constant/g'):
                return "Error: 'constant/g' is not existent while buoyant opiton is selected"
        #paired boundary check: cyclic, porous, etc

    def summarize(self):
        """ summarize the case setup
        """
        print("================= case summary ==================\n")
        print("Solver name: {}".format(self._solverName))
        print("Solver template: {}".format(self._templatePath))
        print("Solver case path: {}".format(self._casePath))
        self._summarizeInternalFields()
        print("Defined mesh boundary names: {}".format(listBoundaryNames(self._casePath)))
        print("Boundary conditions setup:")
        for bc in self._boundarySettings:
            print(bc)
        if self._solverSettings['transient']:
            print("Transient settings:")
            print(self._transientSettings)
        cmdline = self.getSolverCmd()
        # foamJob <solver> & foamLog
        print("Please run the command in new terminal: \n" + cmdline)

    def _summarizeInternalFields(self):
        print("Solver created fields are initialized with value:\n")
        for var in self._solverCreatedVariables:
            f = ParsedParameterFile(self._casePath + os.path.sep + "0" + os.path.sep + var)
            print("    {}:      {}".format(var, f['internalField']))

    def editCase(self):
        "Edit Foam case by platform file explorer like Nautilus on Ubuntu"
        import subprocess
        import sys
        path = self._casePath
        if sys.platform == 'darwin':
            subprocess.Popen(['open', '--', path])
        elif sys.platform.find("linux") == 0:  # python 2: linux2 or linux3; but 'linux' for python3
            subprocess.Popen(['xdg-open', path])
        elif sys.platform == 'win32':
            subprocess.Popen(['explorer', path]) # check_call() will block the python code

    def viewResult(self):
        "view by external program paraview"
        if self._solverSettings['parallel']:
            # if Allrun is excuted, it should reconstruct result
            runFoamApplication(['reconstructPar'], self._casePath)
        runFoamApplication(['paraFoam'], self._casePath)

    def exportResult(self):
        "export to VTK legacy format, ascii or binary format"
        if self._solverSettings['parallel']:
            runFoamApplication(['reconstructPar'],  self._casePath)
        if os.path.exists(self._casePath + os.path.sep + "VTK"):
            shutil.rmtree(self._casePath + os.path.sep + "VTK")
        #pointSetName = 'wholeDomain'
        #createRawFoamFile(self._casePath, 'system', 'topoSetDict', getTopoSetDictTemplate(pointSetName, 'pointSet', boundingBox))
        #runFoamCommand(['topoSet', '-case', self._casePath, '-latestTime'])
        #runFoamCommand(['foamToVTK', '-case', self._casePath, '-latestTime', '-pointSet', pointSetName])
        runFoamApplication(['foamToVTK', '-latestTime'], self._casePath)
        #search for *.vtk
        import glob
        vtk_files = glob.glob(self._casePath + os.path.sep + "VTK" + os.path.sep + "*.vtk")
        if len(vtk_files) >= 1:
            print("only one file name with full path is expected for the result vtk file")
            return vtk_files[-1]


    ###########################################################################
    def getSolverCommand(self):
        if os.path.exists(self._casePath + os.path.sep + "Allrun"):
            return  makeRunCommand("./Allrun", self._casePath)
        else:
            return makeRunCommand(self.getSolverName(), self._casePath)

    def getSolverName(self):
        return _getSolverName(self._solverSettings)

    def getFoamTemplate(self):
        """
        zipped template 'simpleFoam' works with OpenFOAM 3.0+ and 2.x
        In case PyFoam can not load the case file as dict
        """
        solver_name = self._solverName
        if self._solverSettings['turbulenceModel'] in LES_turbulence_models:
            using_LES = True
        else:
            using_LES = False

        script_path = os.path.dirname(os.path.abspath( __file__ ))
        template_path = script_path + os.path.sep + solver_name + "_template_v" + str(getFoamVersion()[0]) + ".zip"
        if not os.path.exists(template_path):
            if solver_name in _RAS_solver_templates:
                template = _RAS_solver_templates[solver_name]
            elif using_LES and solver_name in _LES_turbulenceModel_templates:
                template = _LES_solver_templates[solver_name]
            else:
                raise Exception('No tutorial template is found for solver: ' + solver_name)
            if template:
                return getFoamDir() + os.path.sep + template
            else:
                raise Exception('No tutorial template is registered for solver: ' + solver_name)
        else:
            return template_path # case folder zipped with version number

    def getSolverCreatedVariables(self):
        return set(getVariableList(self._solverSettings))

    def _createInitVarables(self):
        """assuming a empty 0 folder, and all variable file will be created from scratch
        a default fvScheme is needed for each var, but * can be used to match all
        """
        vars = self.getSolverCreatedVariables()
        casePath = self._casePath
        initFolder = casePath + os.path.sep + "0"
        if not os.path.exists(initFolder):
            if _debug:
                print('Warning: init folder (time 0) does not exist, create a new one')
            os.makedirs(initFolder) # mkdir -p

        if _debug: print("Info: initialize solver created fields (variables): ", vars)
        #for f in listVariablesInFolder(initFolder):
        #    #delete not necessary variables
        for v in vars:
            fname = casePath + os.path.sep + "0" + os.path.sep + v
            #clean this file contect if existent
            lines = [
            "dimensions  [0 0 0 0 0 0 0];\n", "\n",
            "internalField uniform 0;\n", "\n",  # will be set again in caseBuilder
            'boundaryField\n', "{\n", "\n", "}\n",
            ]
            if v.split('.')[0] in set([ 'U', 'pointDisplacement']):
                createRawFoamFile(casePath, '0', v, lines, 'volVectorField')
            elif v.split('.')[0] in set(['p', 'p_rgh', 'T', 'alphat', 'alpha',
                                         'k', 'epsilon', 'omega', 'nut', 'nuSgs', 'nuTilda',]):
                createRawFoamFile(casePath, '0', v, lines, 'volScalarField')
            else:
                print("Info: variable {} is left unchanged".format(v))

            f = ParsedParameterFile(fname)
            if v.split('.')[0] == 'U':
                f['dimensions'] = "[0 1 -1 0 0 0 0]"
                f['internalField'] = "uniform (0 0 0)"
            elif v == 'p' or v == "p_rgh":
                if self._solverSettings['compressible']:
                    f['dimensions'] = "[1 -1 -2 0 0 0 0]"  # SI unit [Pa], m-1, kg, s-2
                    f['internalField'] = 'uniform 101000'  # P_rgh for heat transfer P must be > 0
                else:
                    f['dimensions'] = "[0 2 -2 0 0 0 0]"  # it is different from SI unit [Pa], m-1, kg, s-2
                    f['internalField'] = 'uniform 0'  # should be zero for incompressible flow
            elif v == "pointDisplacement":  # dynamicMeshing 6DoF motion solver needs this var
                f['dimensions'] = "[0 1 0 0 0 0 0]"
                f['internalField'] = 'uniform (0 0 0)'
            # thermal variable
            elif v.split('.')[0] == 'T':
                f['dimensions'] = "[0 0 0 1 0 0 0]"
                f['internalField'] = "uniform 300"
            elif v.split('.')[0] == 'alphat':  # thermal turbulence viscosity/diffusivity
                f['dimensions'] = "[1 -1 -1 0 0 0 0]"
                f['internalField'] = 'uniform 0' #
            # default internal fields: see $FOAM_DIR/tutorials/compressible/rhoPimpleFoam/ras/cavity/0
            elif v.split('.')[0] == 'k':
                f['dimensions'] = "[ 0 2 -2 0 0 0 0 ]"
                f['internalField'] = 'uniform 0.00325'
            elif v.split('.')[0] == 'epsilon':
                f['dimensions'] = "[ 0 2 -3 0 0 0 0 ]"
                f['internalField'] = 'uniform 0.000765'
            elif v.split('.')[0] == 'omega':
                f['dimensions'] = "[0 0 -1 0 0 0 0]"  # same unit with Frequency
                f['internalField'] = "uniform 2.6"
            elif v.split('.')[0] == 'nut' or v == 'nuTilda' or v == 'nuSgs':  # mut is needed for openfoam 2.x, foam-ext
                f['dimensions'] = "[0 2 -1 0 0 0 0]"
                f['internalField'] = "uniform 0"
            # variables for multiphase model
            elif v == 'alpha':  # specie fraction
                f['dimensions'] = "[0 0 0 0 0 0 0]"
                f['internalField'] = 'uniform 0'  # non dimensional (0~1)
            else:
                print("variable {} is not recognized and dimension is left unchanged".format(v))

            f.writeFile()

    @property
    def solverSettings(self):
        return self._solverSettings
    def setupSolverSettings(self, settings):
        if settings and isinstance(value, dict):
            self._solverSettings = settings
            self._solverName = getSolverName(self._solverSettings)
            self._solverCreatedVariables = getSolverCreatedVariables(self._solverSettings)

    @property
    def parallelSettings(self):
        return self._parallelSettings
    def setupParallelSettings(self, settings):
        # copy if not existent dic file`decomposeParDict`
        # OpenFOAM 3.0 + can have dict file in other place
        # `decomposeParDict -decomposeParDict dictPath`

        f = self._casePath + os.path.sep + 'decomposeParDict'
        if not os.path.exists(f):
            createRawFoamFile(self._casePath, 'system', 'decomposeParDict',
                              getDecomposeParDictTemplate(4, 'scotch'))
        if settings and isinstance(settings, dict):
            d = ParsedParameterFile(f)
            if 'method' in settings:
                d['method'] = settings['method']
            if 'method' in settings:
                d['numberOfSubdomains'] = settings['numberOfSubdomains']
            d.writeFile()
        #runCommand('decomposePar -case {}'.format(self._casePath)) # maybe too early to call it now

    @property
    def fluidProperties(self):
        return self._fluidProperties
    @fluidProperties.setter
    def fluidProperties(self, value):
        if value and isinstance(value, dict):
            self._fluidProperties = value
        else:
            print("set a invalid fluid property, null or not dict")
        if _debug: print(self._fluidProperties)

    def setupFluidProperties(self):
        """ also control the multiphase properties for incompressible flow
        copy dict file from case template is the preferred way for Nonewtonian fluid
        tutorials/compressible/rhoSimpleFoam: thermophysicalProperties
        tutorials/incompressible/simpleFoam, nonNewtonianIcoFoam: transportProperties
        foam-extend: tutorials/viscoelastic/viscoelasticFluidFoam/Oldroyd-B/constant/viscoelasticProperties
        OpenFoam 3.0, does not need 'nu' before dimension signature, but back-compactible
        """
        self.setupTransportProperties()

    def setupTransportProperties(self):
        case = self._casePath
        solver_settings = self._solverSettings
        assert solver_settings['compressible'] == False

        lines = ['transportModel  Newtonian;\n']
        # PyFOAM 0.6.6 still can not parse `nu    [ 0 2 -1 0 0 0 0 ] 1e-6;` without extra `nu` before dims
        # therefore, can not use `ParsedParameterFile(case + "/constant/transportProperties")`
        if solver_settings['nonNewtonian']:
            print('Warning: nonNewtonian case setup is not implemented, please edit dict file directly')
        else:
            for k in self._fluidProperties:
                if k in set(['nu', 'kinematicViscosity']):
                    viscosity = self._fluidProperties[k]
                    lines.append('nu              nu [ 0 2 -1 0 0 0 0 ] {};\n'.format(viscosity))
                elif k in set(['rho', 'density']):
                    density = self._fluidProperties[k]
                    lines.append('rho              rho [ -3 1 0 0 0 0 0 ] {};\n'.format(density))
                else:
                    print("Warning:unrecoginsed fluid properties: {}".format(k))
            if _debug:
                print("Viscosity settings in constant/transportProperties")
                print(lines)

        createRawFoamFile(case, "constant", "transportProperties", lines)

    @property
    def gravity(self):
        return self._solverSettings['gravity']
    @gravity.setter
    def setGravityProperties(self, value):
        """Assume gravity vector is available in solver settings, creat `constant/g` if not existent
        """
        if value and len(value) == 3:
            self._solverSettings['gravity'] = value

    def setupGravityProperties(self):
        """Assume gravity vector is available in solver settings, creat `constant/g` if not existent
        """
        value = self._solverSettings['gravity']
        lines = ["dimensions    [0 1 -2 0 0 0 0];\n", "value    ({} {} {}); \n".format(value[0], value[1], value[2])]
        createRawFoamFile(self._casePath, "constant", "g", lines, "uniformDimensionedVectorField")

    @property
    def internalFields(self):
        return self._internalFields
    @internalFields.setter
    def internalFields(self, fields=None):
        if fields and isinstance(fields, dict):
            self._internalFields = fields  # mapping type like dict
        else:
            print('Warning: only dict with variable name as key is accepted to init internal fields')

    def setupInternalFields(self, fields=None):
        # support a short name in openfoam dict file like 'T' and long name like 'Temperature'
        if fields and isinstance(fields, dict):
            self.internalFields = fields

        for k in self._internalFields:
            if k.split('.')[0] == 'U': # multiphase has variable name 'U.air'
                var = k
            elif k.split('.')[0] in set(['velocity', 'Velocity']):
                var = k.replace(k.split('.')[0], 'U')
                self._internalFields[var] = self._internalFields[k]
            elif k.lower() in set(['p', 'p_rgh']): # mixture has only one pressure
                var = k
            elif k.lower() in set(['pressure']):
                var = 'p'
                self._internalFields[var] = self._internalFields[k]
            elif k == 'pointDisplacement': # dynamic meshing 6DOF
                var =k
            # turbulence variables
            elif k.split('.')[0] in set(['nut', 'turbulenceviscosity', 'turbulenceViscosity', 'TurbulenceViscosity']):
                var = k.replace(k.split('.')[0], 'nut')
                self._internalFields[var] = self._internalFields[k]
            elif k.split('.')[0] in set(['k', 'turbulenceEnergy', 'TurbulenceEnergy', 'turbulenceKineticEnergy']):
                var = k.replace(k.split('.')[0], 'k')
                self._internalFields[var] = self._internalFields[k]
            elif k.lower() in set(['epsilon', 'turbulencedissipationrate']):
                var = 'epsilon'
            elif k.split('.')[0] == 'omega':
                var = k
            elif k.lower() in set(['omega', 'specificdissipationrate']):
                var = 'omega'
            # thermal flow varaible
            elif k.split('.')[0] in set(['T', 'temperature', 'Temperature']):
                var = k.replace(k.split('.')[0], 'T')
                self._internalFields[var] = self._internalFields[k]
            elif k.lower() in set(['alphat', 'turbulencethermaldiffusivity']):
                var = 'alphat'
            elif k.split('.')[0] == "alphat":
                var = k
            # multiple phase varaible
            elif k.lower().split('.')[0] == 'alpha':
                var = k
            elif k.lower() in set(['alpha', 'phaseratio']):
                var = 'alpha'
            else:
                var = k
            # update the initial internal field
            value = self._internalFields[var]
            if var in self._solverCreatedVariables:
                f = ParsedParameterFile(self._casePath + "/0/" + var)
                f["internalField"] = formatValue(value)
                f.writeFile()
            else:
                print("Warning: variable:{} is not recognised/created thus ignored in setupInternalFields()".format(k))

    """
    see: http://cfd.direct/openfoam/user-guide/controldict/
    startFrom   Controls the start time of the simulation.
    - firstTime: Earliest time step from the set of time directories.
    - startTime: Time specified by the startTime keyword entry.
    - latestTime: Most recent time step from the set of time directories.
    startTime   Start time for the simulation with startFrom startTime;
    stopAt  Controls the end time of the simulation.
    - endTime : Time specified by the endTime keyword entry.
    - writeNow : Stops simulation on completion of current time step and writes data.
    - noWriteNow: Stops simulation on completion of current time step and does not write out data.
    - nextWrite: Stops simulation on completion of next scheduled write time, specified by writeControl.
    """
    @property
    def transientSettings(self):
        return self._transientSettings
    @transientSettings.setter
    def transientSettings(self, transientSettings):
        if transientSettings:
            self._transientSettings = transientSettings

    def setupTransientSettings(self, tSettings):
        if tSettings:
            self.transientSettings = tSettings
        if self._transientSettings:
            f = ParsedParameterFile(self._casePath + "/system/controlDict")
            f["startTime"] = self._transientSettings["startTime"]
            f["endTime"] = self._transientSettings["endTime"]
            f["deltaT"] = self._transientSettings["timeStep"]
            f["writeInterval"] = self._transientSettings["writeInterval"]
            f.writeFile()
        else:
            print("transientSettings is None")

    def setupDynamicMeshingProperties(self):
        """ must relies on template setup,
        also the 'pointDisplacement' boundary need setup
        """
        pass

    ################################## solver control #####################################
    def setupSolverControl(self):
        # pRefValue must be set, since some case does not contains such requested
        # residual control, currently all default from tutorial case setup, but can be setup
        pass

    def setupRelaxationFactors(self, pressure_factor=0.1, velocity_factor=0.1, other_factor=0.3):
        # '.*' apply to all variables
        f = ParsedParameterFile(self._casePath + "/system/fvSolution")
        f["relaxationFactors"]["fields"] = {} # this `fields` may not exist
        f["relaxationFactors"]["fields"]["p"] = pressure_factor
        f["relaxationFactors"]["equations"] = {}
        f["relaxationFactors"]["equations"]["U"] = velocity_factor
        f["relaxationFactors"]["equations"]['"(k|epsilon|omega|nut|f|v2)"'] = other_factor
        f["relaxationFactors"]["equations"]['\".*\"'] = other_factor
        f.writeFile()

    def setupPressureReference(self, pRefValue, pRefCell=0):
        f = ParsedParameterFile(self._casePath + "/system/fvSolution")
        for algo in _supported_algorithms:
            if algo in f:
                f[algo]["pRefValue"] = pRefValue
                f[algo]["pRefCell"] = pRefCell
        f.writeFile()

    def setupNonOrthogonalCorrectors(self, nTimes):
        ## important for netgen mesh with only tetragal cell to convergent
        f = ParsedParameterFile(self._casePath + "/system/fvSolution")
        for algo in _supported_algorithms:
            if algo in f:
                f[algo]['nNonOrthogonalCorrectors'] = nTimes
        f.writeFile()

    def setupResidualControl(self, pResidual, UResidual, other = 0.001):
        f = ParsedParameterFile(self._casePath + "/system/fvSolution")
        for algo in _supported_algorithms:
            if algo in f:
                f[algo]['residualControl'] = {'p': pResidual, 'U': pResidual, '\".*\"': other}
        f.writeFile()

    ############################## non public API ######################################

    def initBoundaryConditions(self):
        """ set all boundary as wall if no boundary setting is specified for this boundary
        write default value in 'case/0/p' 'case/0/U' and turbulence related case/0/var
        in polyMesh/boundary  "defaultFaces" must be wall type,
        but mesh conversion tool does not rename
        """
        bc_names = listBoundaryNames(self._casePath)
        specified_bc_names = set([bc['name'] for bc in self._boundarySettings])
        unspecified_bc_names = set(bc_names) - specified_bc_names
        #for bc in unspecified_bc_names:
        self._initVelocityBoundaryAsWall()
        self._initPressureBoundaryAsWall()
        if 'pointDisplacement' in self._solverCreatedVariables:
            self._initPressure_rghAsWall()
        #thermal boundary init will be called in derived class
        for bc in bc_names:
            _default_wall_dict = {'name': bc, 'type': 'wall', 'subtype': 'fixed'}
            self.setupWallTurbulence(_default_wall_dict, self._turbulenceProperties)

    def _initVelocityBoundaryAsWall(self):
        """incompressible flow has different wall_function:
        """
        f = ParsedParameterFile(self._casePath + "/0/U")
        for bcName in listBoundaryNames(self._casePath):
            f["boundaryField"][bcName]={}
            f["boundaryField"][bcName]["value"]="uniform (0 0 0)"
            f["boundaryField"][bcName]["type"]="fixedValue"
        f.writeFile()

    def _initPressureBoundaryAsWall(self):
        """ shared by compressible flow, porous, nonNewtonian flow
        """
        f = ParsedParameterFile(self._casePath + "/0/p")
        if 'p_rgh' in self._solverCreatedVariables:
            self._initPressure_rghAsWall()
            for bcName in listBoundaryNames(self._casePath):
                f["boundaryField"][bcName]={"type": "calculated", "value": "$internalField"}
                f["boundaryField"][bcName]["type"]="zeroGradient"  # zeroGradient needs no 'value'
        else:
            for bcName in listBoundaryNames(self._casePath):
                f["boundaryField"][bcName]={"type":"zeroGradient" }
        f.writeFile()

    def _initPressure_rghAsWall(self):
        #
        p_rgh = ParsedParameterFile(self._casePath + "/0/p_rgh")
        for bcName in listBoundaryNames(self._casePath):
            p_rgh["boundaryField"][bcName] = {"type": "fixedFluxPressure", 'rho': 'rhok', 'value': "uniform 0"}
        p_rgh.writeFile()

    def _initPointDisplacementAsWall(self):
        #
        pd = ParsedParameterFile(self._casePath + "/0/pointDisplacement")
        for bcName in listBoundaryNames(self._casePath):
            pd["boundaryField"][bcName] = {"type": "fixedValue", 'value': "uniform (0 0 0)"}
        pd.writeFile()

    #####################################################################

    @property
    def boundaryConditions(self):
        return self._boundarySettings
    @boundaryConditions.setter
    def boundaryConditions(self, boundarySettings=None):
        if boundarySettings and isinstance(boundarySettings, list) and len(boundarySettings)>=1:
            self._boundarySettings = boundarySettings

    def setupBoundaryConditions(self, settings=None):
        if settings and isinstance(settings, list) and len(settings)>=1:
            self.boundaryConditions = settings

        self.initBoundaryConditions()
        if not len(self._boundarySettings):
            print("Error: No boundary condition is defined, please check!")
        for bcDict in self._boundarySettings:
            if bcDict['type'] in supported_interface_types:
               self.setupInterfaceBoundary(bcDict['type'], bdDict['name'])
            else:
                assert bcDict['type'] in supported_boundary_types
                if bcDict['type'] == 'inlet':
                    if bcDict['subtype'] == 'totalPressure':
                        self.setupPressureInletBoundary(bcDict)
                    else:  # massflow or uniformVelocity
                        self.setupVelocityInletBoundary(bcDict)
                elif bcDict['type'] == 'outlet':
                    self.setupOutletBoundary(bcDict)
                elif bcDict['type'] == 'freestream':
                    self.setupFreestreamBoundary(bcDict)
                elif bcDict['type'] == 'wall':
                    self.setupWallBoundary(bcDict)
                else:
                    print("Warning: boundary type: {} is not supported yet and ignored!".format(bcDict['type']))
        # extra or special varialbe setup
        if 'p_rgh' in self._solverCreatedVariables:
            self._setupPressure_rgh()
        # pointDisplacement is implemented into setupWallBoundary, setupInletBoundary, etc

    def setupWallBoundary(self, bcDict):
        """
        special wall boundary in OpenFOAM: fixed/noslip, slip, partialSlip, moving, rough
        for heat transfer: see heat transfer builder
        for multiphase flow: alphaContactAngle
        wall turbulence setup has been done in `initTurbulenceBoundaryAsWall()`
        see: http://openfoam.com/documentation/cpp-guide/html/a11608.html
        """
        bcName = bcDict['name']
        wall_type = bcDict['subtype']
        if "value" in bcDict:
            value = bcDict['value']

        f = ParsedParameterFile(self._casePath + "/0/U")
        f["boundaryField"][bcName] = {}
        if wall_type == 'fixed' or wall_type == 'noSlip':  # noSlip equal to fixed wall
            f["boundaryField"][bcName]["type"] = "fixedValue"
            f["boundaryField"][bcName]["value"] = "uniform (0 0 0)"  # initial value
        elif wall_type == 'rough':  # for turbulence flow only
            f["boundaryField"][bcName]["type"] = "roughWall"
            f["boundaryField"][bcName]["U"] = formatValue(value)
            print("Info: wall function for k, epsilon, omege nu may need set for roughwall")
        elif wall_type == 'slip':  # 100% slipping between solid and fluid
            f["boundaryField"][bcName]["type"] = "slip"
        elif wall_type == 'partialSlip':  # for multiphase flow
            f["boundaryField"][bcName]["type"] = "partialSlip"
            f["boundaryField"][bcName]["valueFraction"] = value
            f["boundaryField"][bcName]["value"] = "uniform (0 0 0)" #
        elif wall_type == "moving":  # translatingWallVelocity  movingWallVelocity
            f["boundaryField"][bcName]["type"] = "movingWallVelocity"
            f["boundaryField"][bcName]["U"] = formatValue(value)
            if self._solverSettings['dynamicMeshing'] and "pointDisplacement" in self._solverCreatedVariables():
                df = ParsedParameterFile(self._casePath + "/0/pointDisplacement")
                df["boundaryField"][bcName] = {'type': "calculated", 'value': 'uniform (0 0 0)'}
                df.writeFile()
        else:
            print("wall boundary: {} is not supported yet".format(wall_type))
        f.writeFile()

        pf = ParsedParameterFile(self._casePath + "/0/p")
        pf["boundaryField"][bcName] = {'type': "zeroGradient"}
        pf.writeFile()

        if 'turbulenceSettings' in bcDict:
            turbulenceSettings = bcDict['turbulenceSettings']
        else:
            turbulenceSettings = self._turbulenceProperties
        self.setupWallTurbulence(bcDict, turbulenceSettings)

    def setupInterfaceBoundary(self, bcDict):
        # freestream is kind of like interface boundary type like wall, always uniform $internalField
        ## empty:  2D as a special case of 3D mesh with one layer of mesh with front and back as empty
        # symmetry: for any (non-planar) patch which uses the symmetry plane (slip) condition.
        # "cyclicAMI" for Rotating, "cyclicAMI" slidingMesh
        # "processor": for MPI parallel calculation, user needs not setup
        case = self._casePath
        interface_type = bcDict['subtype']
        boundary_name = bcDict['name']
        var_list = self._solverCreatedVariables
        for var in var_list:
            f = ParsedParameterFile(case + "/0/" + var)
            if interface_type == "empty" or interface_type == "2Dinerface":  # 2D case single layer extruded 3D emsh
                f["boundaryField"]["frontAndBack"] = {}
                f["boundaryField"]["frontAndBack"]["type"] = "empty"
                # axis-sym 2D case, axis line is also empty type
            elif interface_type == "symmetryPlane" or interface_type == "symmetry":
                f["boundaryField"][boundary_name] = {}
                f["boundaryField"][boundary_name]["type"] = interface_type
            elif interface_type == "cyclic": # also named as `periodic` in ansys Fluent
                f["boundaryField"][boundary_name] = {}
                f["boundaryField"][boundary_name]["type"] = "cyclic"
                # todo: chack pairing of wedge and cyclic boundary
            elif interface_type == "wedge":
                f["boundaryField"][boundary_name] = {}
                f["boundaryField"][boundary_name]["type"] = "wedge"  # axis-sym
            elif interface_type == "coupled":
                raise Exception('Boundary or patch type {} is not supported'.format(interface_type))
                # todo: need to pair with other analysis
            else:
                raise Exception('Boundary or patch type {} is not supported'.format(interface_type))
            f.writeFile()

    def pairCyclicBoundary(self, type, names):
        """
        http://www.cfdsupport.com/OpenFOAM-Training-by-CFD-Support/node108.html
        inGroups 1(cyclicAMI); neighbourPatch <ref to the paired patch name>
        #boundary file needs to be modified,
        'transform' = roational, rotationAxis. rotatingCentre
        'transform' = translational; separationVector =
        """
        pass

    def _setupPressure_rgh(self):
        # Pseudo hydrostatic pressure [Pa] : p - rho * g * height
        # only for buoyant flow in heat transfer in compressible and incompressible flow
        # /opt/openfoam4/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.orig/porous
        # /opt/openfoam4/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.orig
        # prghTotalPressure,
        # buoyantPressure for heat exchanger wall, value uniform 0;
        # source doc: http://www.openfoam.com/documentation/cpp-guide/html/a02111.html#details

        f = ParsedParameterFile(self._casePath + "/0/p_rgh")
        for bcDict in self._boundarySettings:
            bc = bcDict['name']
            subtype = bcDict['subtype']
            if bcDict['type'] in set(['inlet', 'outlet']):
                if bcDict['subtype'] in set(["staticPressure", "pressure"]):
                    f["boundaryField"][bc] = {'type': 'prghPressure', 'p': formatValue(bcDict['value'])}
                if bcDict['subtype'] == "totalPressure":
                    f["boundaryField"][bc] = {'type': 'prghTotalPressure', 'p0': formatValue(bcDict['value'])}
                else:
                    f["boundaryField"][bc] = {'type': 'fixedFluxPressure', 'gradient':'uniform 0', 'value': "$internalField"}
            elif bcDict['type'] == 'freestream':
                f["boundaryField"][bc] = {'type': 'fixedValue', 'value': "$internalField"} #todo: check freestreamPressure
            elif bcDict['type'] == 'interface':
                f["boundaryField"][bc] = {'type': subtype}
            elif bcDict['type'] == 'wall':
                f["boundaryField"][bc] = {'type': 'fixedFluxPressure', 'value': 'uniform 0'}
            else:
                print("Warning: boundary type: {} is not supported in p_rgh thus ignored!".format(bcDict['type']))
        f.writeFile()

    def setupPressureInletBoundary(self, bcDict):
        # value is MPa in FreeCAD, but Pa is needed in OpenFOAM
        bcName = bcDict['name']
        inlet_type = bcDict['subtype']
        value = bcDict['value']

        pf = ParsedParameterFile(self._casePath + "/0/p")
        pf["boundaryField"][bcName] = {}
        Uf = ParsedParameterFile(self._casePath + "/0/U")
        Uf["boundaryField"][bcName] = {}

        if 'p_rgh' in self._solverCreatedVariables:
            pf["boundaryField"][bcName] = {'type': 'calculated', 'value': "$internalField"}
        else:
            if inlet_type == "totalPressure":
                pf["boundaryField"][bcName]["type"] = 'totalPressure'
                pf["boundaryField"][bcName]["p0"] = 'uniform {}'.format(value)
                #pf["boundaryField"][bcName]["gamma"] = 0  # for air: 1 .4
                pf["boundaryField"][bcName]["value"] = "$internalField"  # initial value
            else:  #
                pf["boundaryField"][bcName]["type"] = 'fixedValue'
                pf["boundaryField"][bcName]["value"] = "uniform {}".format(value)
        pf.writeFile()

        # velocity intial value is default to wall, uniform 0, so it needs to change
        Uf["boundaryField"][bcName]["type"] = "pressureInletOutletVelocity"
        Uf["boundaryField"][bcName]["value"] ="uniform (0 0 0)"  # initial value
        Uf.writeFile()

        if 'turbulenceSettings' in bcDict:
            turbulenceSettings = bcDict['turbulenceSettings']
        else:
            turbulenceSettings = self._turbulenceProperties
        self.setupInletTurbulence(bcDict, turbulenceSettings)

    def setupVelocityInletBoundary(self, bcDict):
        """ direction: by default, normal to inlet boundary
        """
        bcName = bcDict['name']
        inlet_type = bcDict['subtype']
        value = bcDict['value']

        # velocity intial value is default to wall: uniform (0,0,0)
        Uf = ParsedParameterFile(self._casePath + "/0/U")
        Uf["boundaryField"][bcName] = {}
        if inlet_type == "massFlowRate":  # compressible flow only?
            Uf["boundaryField"][bcName]["type"] = "flowRateInletVelocity"
            Uf["boundaryField"][bcName]["massFlowRate"] = value  # kg/s
            Uf["boundaryField"][bcName]["rho"] = "rho"
            Uf["boundaryField"][bcName]["rhoInlet"] = "rho"
            Uf["boundaryField"][bcName]["value"] = "$internalField"
        elif inlet_type == "volumetricFlowRate":
            Uf["boundaryField"][bcName]["type"] = "flowRateInletVelocity"
            Uf["boundaryField"][bcName]["volumetricFlowRate"] = value  # m3/s
            Uf["boundaryField"][bcName]["value"] = "uniform (0 0 0)"
        elif inlet_type == "uniformVelocity":
            assert len(value) == 3  # velocity must be a tuple or list with 3 components
            Uf["boundaryField"][bcName]["type"] = "fixedValue"
            Uf["boundaryField"][bcName]["value"] = formatValue(value)
        else:
            print(inlet_type + " is not supported as inlet boundary type")
        Uf.writeFile()

        pf = ParsedParameterFile(self._casePath + "/0/p")
        if 'p_rgh' in self._solverCreatedVariables:
            pf["boundaryField"][bcName] = {'type': 'calculated', 'value': "$internalField"}
        else:
            pf["boundaryField"][bcName] = {}
            pf["boundaryField"][bcName]["type"] = "zeroGradient"
            pf["boundaryField"][bcName]["value"] = "uniform 0"
        pf.writeFile()
        #
        if 'turbulenceSettings' in bcDict:
            turbulenceSettings = bcDict['turbulenceSettings']
        else:
            turbulenceSettings = self._turbulenceProperties
        self.setupInletTurbulence(bcDict, turbulenceSettings)

    def setupOutletBoundary(self, bcDict):
        """supported_outlet_types: outFlow, pressureOutlet
        pressureOutlet, for exit to background static pressure
        outFlow: corresponding to velocityInlet/FlowRateInlet
        self.supported_outlet_types = set([])
        http://cfd.direct/openfoam/user-guide/boundaries/
        inletOutlet: Switches U and p between fixedValue and zeroGradient depending on direction of U
        pressureInletOutletVelocity: Combination of pressureInletVelocity and inletOutlet
        """
        bcName = bcDict['name']
        outlet_type = bcDict['subtype']
        value = bcDict['value']

        pf = ParsedParameterFile(self._casePath + "/0/p")
        pf["boundaryField"][bcName] = {}
        if outlet_type == "totalPressure":
            pf["boundaryField"][bcName]["type"] = 'totalPressure'
            pf["boundaryField"][bcName]["p0"] = 'uniform {}'.format(value)
            pf["boundaryField"][bcName]["gamma"] = 0  # what?  1 .4
            pf["boundaryField"][bcName]["value"] = "$internalField"  # initial value
        elif outlet_type == "staticPressure":
            pf["boundaryField"][bcName]["type"] = "fixedValue"  # totalPressure
            pf["boundaryField"][bcName]["value"] = "uniform {}".format(value)
        elif outlet_type == "massFlowRate" or outlet_type == "volumetricFlowRate":
            #pf["boundaryField"][bcName]["type"] = "zeroGradient"
            print("Error: massFlowRate  and volumetricFlowRate not yet supported")
            return
        elif outlet_type == "outFlow": #PressureInletVelocityOutlet
            pf["boundaryField"][bcName]["type"] = "outletInlet"
            #pf["boundaryField"][bcName]["outletValue"] = "uniform 0"  # not sure! check
            pf["boundaryField"][bcName]["value"] ="$internalField"
        else:
            #default to zeroGradient, for velocityOutlet
            pf["boundaryField"][bcName]["type"] = "zeroGradient"
            print("pressure bundary default to zeroGradient for outlet type '{}' ".format(outlet_type))
        pf.writeFile()
        # velocity intial value is default to wall, uniform 0, so it needs to change
        Uf = ParsedParameterFile(self._casePath + "/0/U")
        Uf["boundaryField"][bcName] = {}
        if outlet_type == "totalPressure" or outlet_type == "staticPressure" :
            Uf["boundaryField"][bcName]["type"] = "pressureInletOutletVelocity"  #
            Uf["boundaryField"][bcName]["value"] ="uniform (0 0 0)"  # set as initial value only
        elif outlet_type == "uniformVelocity":
            Uf["boundaryField"][bcName]["type"] = "fixedValue"
            Uf["boundaryField"][bcName]["value"] = formatValue(value)
        elif outlet_type == "outFlow":
            Uf["boundaryField"][bcName]["type"] = "inletOutlet"
            Uf["boundaryField"][bcName]["outletValue"] ="uniform (0 0 0)"
            Uf["boundaryField"][bcName]["value"] ="$internalField"
            #Uf["boundaryField"][bcName]["inletValue"] ="$internalField"
        else:  # "massFlowRate" "volumetricFlowRate"
            print("velocity bundary set to inletOutlet for outlet type '{}' ".format(outlet_type))
            Uf["boundaryField"][bcName]["type"] = "zeroGradient"
        Uf.writeFile()
        #
        if 'turbulenceSettings' in bcDict:
            turbulenceSettings = bcDict['turbulenceSettings']
        else:
            turbulenceSettings = self._turbulenceProperties
        self.setupOutletTurbulence(bcDict, turbulenceSettings)


    def setupFreestreamBoundary(self, bcDict):
        """ freestream is the known velocity in far field of open space
        see: tutorials/incompressible/simpleFoam/airFoil2D/0.org/
        """
        print('internalField must be set for each var in 0/ folder')
        bcName = bcDict['name']
        value = bcDict['value']
        
        if bcDict['type'] == "freestreamPressure" or bcDict['type'] == "freestream":
            f = ParsedParameterFile(self._casePath + "/0/p")
            f["boundaryField"][bcName] = {}
            f["boundaryField"][bcName]["type"] = "freestreamPressure"
            f.writeFile()
            #
            f = ParsedParameterFile(self._casePath + "/0/U")
            f["boundaryField"][bcName] = {}
            f["boundaryField"][bcName]["type"] = "freestream"
            f["boundaryField"][bcName]["value"] = formatValue(value)
            f.writeFile()
        elif bcDict['type'] == "freestreamVelocity":
            f = ParsedParameterFile(self._casePath + "/0/p")
            f["boundaryField"][bcName] = {}
            f["boundaryField"][bcName]["type"] = "freestream"
            f.writeFile()
            #
            f = ParsedParameterFile(self._casePath + "/0/U")
            f["boundaryField"][bcName] = {}
            f["boundaryField"][bcName]["type"] = "freestreamVelocity"
            f["boundaryField"][bcName]["value"] = formatValue(value)
            f.writeFile()
        else:
            print( bcDict['type'] + "is not yet supported")

        # todo: freestream turbulence may need extra setting up
        if 'turbulenceSettings' in bcDict:
            turbulenceSettings = bcDict['turbulenceSettings']
        else:
            turbulenceSettings = self._turbulenceProperties
        self.setupFreestreamTurbulence(bcDict, turbulenceSettings)


    ###################################################################################

    @property
    def turbulenceProperties(self):
        return self._turbulenceProperties

    @turbulenceProperties.setter
    def turbulenceProperties(self, tSettings=None):
        """ see list of turbulence model: http://www.openfoam.org/features/turbulence.php
        OpenFoam V3.0 has unified turbuence setup dic, incompatible with 2.x
        currently only some common RAS models are settable by this script
        """
        if tSettings:
            self._turbulenceProperties = tSettings
            if _debug:
                print(self._turbulenceProperties)

    def setupTurbulenceProperties(self, turbulenceProperties=None):
        """ see list of turbulence model: http://www.openfoam.org/features/turbulence.php
        OpenFoam V3.0 has unified turbuence setup dic, incompatible with 2.x
        currently only some common RAS models are settable by this script
        """
        if turbulenceProperties:
            self.turbulenceProperties = turbulenceProperties
        case = self._casePath
        turbulanceModelName = self._turbulenceProperties['name']
        fname = case + "/constant/turbulenceProperties"

        if not os.path.exists(fname):
            lines = ["simulationType     laminar;"]
            createRawFoamFile(case, "constant", "turbulenceProperties", lines)
        #PyFoam 0.6.6 can not parse some dict for OpenFOAM 3.0+
        f = ParsedParameterFile(case + "/constant/turbulenceProperties")

        if turbulanceModelName in set(['laminar']):
            f['simulationType'] = "laminar"
            if 'RAS' in f:
                del f['RAS']  # clear this content
        elif turbulanceModelName in RAS_turbulence_models:
            f['simulationType'] = "RAS"
            if getFoamVariant() == "OpenFOAM" and getFoamVersion()[0] >= 3:
                f['RAS'] = {'RASModel': turbulanceModelName, 'turbulence': "on", 'printCoeffs': "on"}
                if turbulanceModelName in kEpsilon_models and 'kEpsilonCoeffs' in self._turbulenceProperties:
                    f['kEpsilonCoeffs'] = self._turbulenceProperties['kEpsilonCoeffs']
                # need check again, other model may need diff parameter
            else:
                fRAS = ParsedParameterFile(case + "/constant/RASProperties")
                fRAS['RASModel'] = turbulanceModelName
                fRAS['turbulence'] = "on"
                fRAS['printCoeffs'] = "on"
                # modify coeff is it is diff from default (coded into source code)
                if turbulanceModelName == "kEpsilon" and 'kEpsilonCoeffs' in self._turbulenceProperties:
                    fRAS['kEpsilonCoeffs'] = self._turbulenceProperties['kEpsilonCoeffs']
                fRAS.writeFile()

        elif turbulanceModelName in LES_turbulence_models:
            if getFoamVariant() == "OpenFOAM" and getFoamVersion()[0] >= 3:
                # all LES model setup is done in file 'turbulenceProperties'
                if not os.path.exists(fname):
                    createRawFoamFile(case, "constant", "turbulenceProperties",
                                        getLeSTurbulencePropertiesTemplate(turbulanceModelName))
            else:
                print("Assumeing LESTurbulenceProperties is copied from template file for OpenFOAM 2.x or foam-extend")
        else:
            print("Turbulence model {} is not recognised or implemented".format(turbulanceModelName))
        f.writeFile()

    def listTurbulenceVarables(self):
        """specific turbulenceModel name:  kEpsilon, kOmegaSST, SpalartAllmaras
        SpalartAllmaras: incompressible/simpleFoam/airFoil2D
        'mut', used in compressible fluid
        """
        solverSettings = self._solverSettings
        return getTurbulenceVariables(solverSettings)

    def setupWallTurbulence(self, bcDict, turbulenceSettings):
        """ set all boundary condition to default: wall, also set wallFunction for each var
        for y+<1 (viscous sublayer) turbulent field = 0, if y+>30, wall function should be used
        - compressible flow has suffex for nut: see tutorial of `rhoSimpleFoam`
        - diff turbulence modle has diff wall function: nutSpalartAllmarasWallFunction
        - rough wall has speicial wall function: not supported yet
        - wall function can be related with velocity U, not supported yet
        - low Reynolds model has suffex `lowRe`: supported
        #includeEtc "caseDicts/setConstraintTypes"
        compressible cases have different wall functions with a suffix for OpenFOAM 2.x
        """
        case = self._casePath
        turbulence_var_list = self.listTurbulenceVarables()
        if self._solverSettings['compressible']:
            suffix = 'compressible::'  # only `alphat` is needed for OpenFOAM 3.0+
        else:
            suffix = ''
        #print(turbulence_model)
        if self._solverSettings['turbulenceModel'] in lowRe_models:
            print('Warning: low Reynolds wall function is not supported yet')
            lowRe = "LowRe"
            kWallFunction = 'kLowReWallFunction'
        else:
            lowRe = ""
            kWallFunction = 'kqRWallFunction'

        for var in turbulence_var_list:
            f = ParsedParameterFile(case + "/0/" + var)
            # kOmega has nonzero internalField for k, omega and epsilon, set default in CreateInitVarables()
            bcName = bcDict['name']
            # if boundaryType == 'wall' and 'type == 'rough':
            #   print('rough wall boundary is not support yet')
            f["boundaryField"][bcName] = {}
            if var == 'k' or var.find("k.") == 0: #begin with
                f["boundaryField"][bcName]["value"]="$internalField"
                f["boundaryField"][bcName]["type"]=kWallFunction  # depend on turbulence_model!!!

            elif var == 'epsilon' or var.find("epsilon") == 0:
                f["boundaryField"][bcName]["type"]= var + lowRe + "WallFunction"
                f["boundaryField"][bcName]["value"]="$internalField" #

            elif var == 'omega' or var.find("omega") == 0:
                f["boundaryField"][bcName]["type"]= var + lowRe + "WallFunction"
                f["boundaryField"][bcName]["value"]="$internalField"
            # alphaT may be init/setup again in derived class thermal builder
            elif var == "alphat" or var.find("alphat") == 0:
                f["boundaryField"][bcName]["type"] = "calculated"
                f["boundaryField"][bcName]["value"] = "$internalField"

            elif  var == 'nut' or var.find("nut.") == 0:  # nut nutk
                if  turbulenceSettings['name'] in spalartAllmaras_models:
                    f["boundaryField"][bcName]["value"]="uniform 0"
                    f["boundaryField"][bcName]["type"]="nutSpalartAllmarasWallFunction" # or 'nutUSpaldingWallFunction'
                    print('Info: nutSpalartAllmarasWallFunction/nutUSpaldingWallFunction can be used for spalartAllmaras turbulence model')
                else:
                    f["boundaryField"][bcName]["value"]="$internalField"
                    f["boundaryField"][bcName]["type"]="nutkWallFunction"  # calculated for other type

            elif var == "nuTilda" or var.find("nuTilda.") == 0:   # for LES model and SpalartAllmaras RAS models
                f["boundaryField"][bcName]["type"]="zeroGradient"  # 'value' not needed, same for LES
            else:
                print("Warning: turbulent var {} is not recognised thus ignored".format(var))
            f.writeFile()

    def setupInletTurbulence(self, bcDict, turbulenceSettings):
        """ modeled from case tutorials/incompressible/simpleFoam/pipeCyclic/0.org
        available turbulentce spec:
        - Set k and epislon explicitly.
        - Set turbulence intensity and turbulence length scale.
        - Set turbulence intensity and turbulent viscosity ratio
        - Set turbulence intensity and hydraulic diameter
        """
        case = self._casePath
        bcName = bcDict['name']
        turbulence_var_list = self.listTurbulenceVarables()

        if "turbulentIntensity" in turbulenceSettings:
            turbulentIntensity = turbulenceSettings["intensityValue"]
        else:
            turbulentIntensity = 0.05  # 5% default, a reasonable guess
        if "hydrauicDiameter" in turbulenceSettings:
            turbulentMixingLength = 0.5 * turbulenceSettings["lengthValue"]
        else:
            turbulentMixingLength = 0.1  # in metre, half inlet diam/width
        #print(turbulence_var_list)
        for var in turbulence_var_list:
            f = ParsedParameterFile(case + "/0/" + var)
            f["boundaryField"][bcName] = {}
            if var == 'k' or var.find("k.") == 0: #begin with
                f["boundaryField"][bcName]["type"] = "turbulentIntensityKineticEnergyInlet"
                f["boundaryField"][bcName]["intensity"] = turbulentIntensity
                f["boundaryField"][bcName]["value"] = "$internalField"

            elif var == 'epsilon' or var.find("epsilon") == 0:
                f["boundaryField"][bcName]["type"] = "turbulentMixingLengthDissipationRateInlet"
                f["boundaryField"][bcName]["mixingLength"] = turbulentMixingLength
                f["boundaryField"][bcName]["value"] = "$internalField"

            elif var == 'omega' or var.find("omega") == 0:
                # /opt/openfoam4/etc/templates/axisymmetricJet
                f["boundaryField"][bcName]["type"] = "turbulentMixingLengthFrequencyInlet"
                f["boundaryField"][bcName]["mixingLength"] = turbulentMixingLength
                f["boundaryField"][bcName]["value"] = "$internalField"

            elif var == "alphat" or var.find("alphat") == 0:
                f["boundaryField"][bcName]["type"] = {"type": "calculated",  "value": "$internalField"}

            elif  var == 'nut' or var.find("nut.") == 0:
                f["boundaryField"][bcName] = {"type": "calculated",  "value": "$internalField"}

            elif var == "nuTilda" or var.find("nuTilda.") == 0:
                f["boundaryField"][bcName] = {"type": "zeroGradient"} #zeroGradient; for wall, inlet and outlet
            else:
                print("Warning: turbulent var {} is not recognised thus ignored".format(var))
            f.writeFile()

    def setupOutletTurbulence(self, bcDict, turbulenceSettings):
        case = self._casePath
        bcName = bcDict['name']
        turbulence_var_list = self.listTurbulenceVarables()
        for var in turbulence_var_list:
            f = ParsedParameterFile(self._casePath + "/0/" + var)
            f["boundaryField"][bcName] = {}
            if var == 'k' or var.find("k.") == 0: #begin with
                f["boundaryField"][bcName]["type"] = "inletOutlet"
                f["boundaryField"][bcName]["inletValue"] = "$internalField"
                f["boundaryField"][bcName]["value"] = "$internalField"

            elif var == 'epsilon' or var.find("epsilon") == 0:
                f["boundaryField"][bcName]["type"] = "inletOutlet"
                f["boundaryField"][bcName]["inletValue"] = "$internalField"
                f["boundaryField"][bcName]["value"] = "$internalField"

            elif var == 'omega' or var.find("omega") == 0:
                f["boundaryField"][bcName]["type"] = "inletOutlet"
                f["boundaryField"][bcName]["inletValue"] = "$internalField"
                f["boundaryField"][bcName]["value"] = "$internalField"

            elif var == "alphat" or var.find("alphat") == 0:
                f["boundaryField"][bcName]["type"] = "calculated"
                f["boundaryField"][bcName]["value"] = "$internalField"

            elif  var == 'nut' or var.find("nut.") == 0:
                f["boundaryField"][bcName]["type"] = "calculated"
                f["boundaryField"][bcName]["value"] = "$internalField"

            elif var == "nuTilda" or var.find("nuTilda.") == 0:
                f["boundaryField"][bcName]["type"] = "zeroGradient"
                #zeroGradient; same for wall, inlet and outlet
            else:
                print("Warning: turbulent var {} is not recognised thus ignored".format(var))
            f.writeFile()

    def setupFreestreamTurbulence(self, bcDict, turbulenceSettings):
        # see: tutorials/incompressible/simpleFoam/airFoil2D/0.org/  SA model
        # tutorial for other RAS model?  k, omege, epsilon, nut, alphat,  "calculated"?
        case = self._casePath
        bcName = bcDict['name']
        bType = bcDict['type']
        subtype = bcDict['subtype']
        turbulence_var_list = self.listTurbulenceVarables()
        if not "freestreamValue" in turbulenceSettings:
            turbulenceSettings["freestreamValue"] = '$internalField'
        for v in turbulence_var_list:
            f = ParsedParameterFile(self._casePath + "/0/" + v)
            f["boundaryField"][bcName] = {}
            if v.split('.')[0] in set(['nut', 'nuTilda', 'alphat']):
                f["boundaryField"][bcName]["type"] = "freestream"
                f["boundaryField"][bcName]["freestreamValue"] = "uniform {}".format(turbulenceSettings["freestreamValue"])
            elif v.split('.')[0] in set(['k', 'epsilon', 'omega']):
                f["boundaryField"][bcName]["type"] = "calculated"
                f["boundaryField"][bcName]["value"] = "uniform {}".format(turbulenceSettings["freestreamValue"])
            else:
                print("Warning: turbulent var {} is not recognised thus ignored".format(var))
            f.writeFile()

    def setupInterfaceTurbulence(self, bcDict, turbulenceSettings):
        #todo: need double check with tutorials
        bcName = bcDict['name']
        subtype = bcDict['subtype']
        for v in turbulence_var_list:
            f = ParsedParameterFile(self._casePath + "/0/" + v)
            f["boundaryField"][bcName] = {}
            f["boundaryField"][bcName]["type"] = subtype
            f.writeFile()